home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / slatec.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  7.5 KB  |  206 lines
  1. /* rng/slatec.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /**
  21.  
  22. * ======================================================================
  23. * NIST Guide to Available Math Software.
  24. * Source for module RAND from package CMLIB.
  25. * Retrieved from TIBER on Fri Oct 11 11:43:42 1996.
  26. * ======================================================================
  27.       FUNCTION RAND(R)
  28. C***BEGIN PROLOGUE  RAND
  29. C***DATE WRITTEN   770401   (YYMMDD)
  30. C***REVISION DATE  820801   (YYMMDD)
  31. C***CATEGORY NO.  L6A21
  32. C***KEYWORDS  RANDOM NUMBER,SPECIAL FUNCTION,UNIFORM
  33. C***AUTHOR  FULLERTON, W., (LANL)
  34. C***PURPOSE  Generates a uniformly distributed random number.
  35. C***DESCRIPTION
  36. C
  37. C      This pseudo-random number generator is portable among a wide
  38. C variety of computers.  RAND(R) undoubtedly is not as good as many
  39. C readily available installation dependent versions, and so this
  40. C routine is not recommended for widespread usage.  Its redeeming
  41. C feature is that the exact same random numbers (to within final round-
  42. C off error) can be generated from machine to machine.  Thus, programs
  43. C that make use of random numbers can be easily transported to and
  44. C checked in a new environment.
  45. C      The random numbers are generated by the linear congruential
  46. C method described, e.g., by Knuth in Seminumerical Methods (p.9),
  47. C Addison-Wesley, 1969.  Given the I-th number of a pseudo-random
  48. C sequence, the I+1 -st number is generated from
  49. C             X(I+1) = (A*X(I) + C) MOD M,
  50. C where here M = 2**22 = 4194304, C = 1731 and several suitable values
  51. C of the multiplier A are discussed below.  Both the multiplier A and
  52. C random number X are represented in double precision as two 11-bit
  53. C words.  The constants are chosen so that the period is the maximum
  54. C possible, 4194304.
  55. C      In order that the same numbers be generated from machine to
  56. C machine, it is necessary that 23-bit integers be reducible modulo
  57. C 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
  58. C integers be multiplied exactly.  Furthermore, if the restart option
  59. C is used (where R is between 0 and 1), then the product R*2**22 =
  60. C R*4194304 must be correct to the nearest integer.
  61. C      The first four random numbers should be .0004127026,
  62. C .6750836372, .1614754200, and .9086198807.  The tenth random number
  63. C is .5527787209, and the hundredth is .3600893021 .  The thousandth
  64. C number should be .2176990509 .
  65. C      In order to generate several effectively independent sequences
  66. C with the same generator, it is necessary to know the random number
  67. C for several widely spaced calls.  The I-th random number times 2**22,
  68. C where I=K*P/8 and P is the period of the sequence (P = 2**22), is
  69. C still of the form L*P/8.  In particular we find the I-th random
  70. C number multiplied by 2**22 is given by
  71. C I   =  0  1*P/8  2*P/8  3*P/8  4*P/8  5*P/8  6*P/8  7*P/8  8*P/8
  72. C RAND=  0  5*P/8  2*P/8  7*P/8  4*P/8  1*P/8  6*P/8  3*P/8  0
  73. C Thus the 4*P/8 = 2097152 random number is 2097152/2**22.
  74. C      Several multipliers have been subjected to the spectral test
  75. C (see Knuth, p. 82).  Four suitable multipliers roughly in order of
  76. C goodness according to the spectral test are
  77. C    3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
  78. C    2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
  79. C    3146245 = 1536*2048 +  517 = 2**21 + 2**20 + 2**9 + 5
  80. C    2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
  81. C
  82. C      In the table below LOG10(NU(I)) gives roughly the number of
  83. C random decimal digits in the random numbers considered I at a time.
  84. C C is the primary measure of goodness.  In both cases bigger is better.
  85. C
  86. C                   LOG10 NU(I)              C(I)
  87. C       A       I=2  I=3  I=4  I=5    I=2  I=3  I=4  I=5
  88. C
  89. C    3146757    3.3  2.0  1.6  1.3    3.1  1.3  4.6  2.6
  90. C    2098181    3.3  2.0  1.6  1.2    3.2  1.3  4.6  1.7
  91. C    3146245    3.3  2.2  1.5  1.1    3.2  4.2  1.1  0.4
  92. C    2776669    3.3  2.1  1.6  1.3    2.5  2.0  1.9  2.6
  93. C   Best
  94. C    Possible   3.3  2.3  1.7  1.4    3.6  5.9  9.7  14.9
  95. C
  96. C             Input Argument --
  97. C R      If R=0., the next random number of the sequence is generated.
  98. C        If R .LT. 0., the last generated number will be returned for
  99. C          possible use in a restart procedure.
  100. C        If R .GT. 0., the sequence of random numbers will start with
  101. C          the seed R mod 1.  This seed is also returned as the value of
  102. C          RAND provided the arithmetic is done exactly.
  103. C
  104. C             Output Value --
  105. C RAND   a pseudo-random number between 0. and 1.
  106. C***REFERENCES  (NONE)
  107. C***ROUTINES CALLED  (NONE)
  108. C***END PROLOGUE  RAND
  109.       DATA IA1, IA0, IA1MA0 /1536, 1029, 507/
  110.       DATA IC /1731/
  111.       DATA IX1, IX0 /0, 0/
  112. C***FIRST EXECUTABLE STATEMENT  RAND
  113.       IF (R.LT.0.) GO TO 10
  114.       IF (R.GT.0.) GO TO 20
  115. C
  116. C           A*X = 2**22*IA1*IX1 + 2**11*(IA1*IX1 + (IA1-IA0)*(IX0-IX1)
  117. C                   + IA0*IX0) + IA0*IX0
  118. C
  119.       IY0 = IA0*IX0
  120.       IY1 = IA1*IX1 + IA1MA0*(IX0-IX1) + IY0
  121.       IY0 = IY0 + IC
  122.       IX0 = MOD (IY0, 2048)
  123.       IY1 = IY1 + (IY0-IX0)/2048
  124.       IX1 = MOD (IY1, 2048)
  125. C
  126.  10   RAND = IX1*2048 + IX0
  127.       RAND = RAND / 4194304.
  128.       RETURN
  129. C
  130.  20   IX1 = AMOD(R,1.)*4194304. + 0.5
  131.       IX0 = MOD (IX1, 2048)
  132.       IX1 = (IX1-IX0)/2048
  133.       GO TO 10
  134. C
  135.       END
  136.  
  137.   **/
  138.  
  139. #include <config.h>
  140. #include <stdlib.h>
  141. #include <gsl/gsl_rng.h>
  142.  
  143. static inline unsigned long int slatec_get (void *vstate);
  144. static double slatec_get_double (void *vstate);
  145. static void slatec_set (void *state, unsigned long int s);
  146.  
  147. typedef struct
  148.   {
  149.     long int x0, x1;
  150.   }
  151. slatec_state_t;
  152.  
  153. static const long P = 4194304;
  154. static const long a1 = 1536;
  155. static const long a0 = 1029;
  156. static const long a1ma0 = 507;
  157. static const long c = 1731;
  158.  
  159. static inline unsigned long int
  160. slatec_get (void *vstate)
  161. {
  162.   long y0, y1;
  163.   slatec_state_t *state = (slatec_state_t *) vstate;
  164.  
  165.   y0 = a0 * state->x0;
  166.   y1 = a1 * state->x1 + a1ma0 * (state->x0 - state->x1) + y0;
  167.   y0 = y0 + c;
  168.   state->x0 = y0 % 2048;
  169.   y1 = y1 + (y0 - state->x0) / 2048;
  170.   state->x1 = y1 % 2048;
  171.  
  172.   return state->x1 * 2048 + state->x0;
  173. }
  174.  
  175. static double 
  176. slatec_get_double (void *vstate)
  177. {
  178.   return slatec_get (vstate) / 4194304.0 ;
  179. }
  180.  
  181. static void
  182. slatec_set (void *vstate, unsigned long int s)
  183. {
  184.   slatec_state_t *state = (slatec_state_t *) vstate;
  185.  
  186.   /* Only eight seeds are permitted.  This is pretty limiting, but
  187.      at least we are guaranteed that the eight sequences are different */
  188.  
  189.   s = s % 8;
  190.   s *= P / 8;
  191.  
  192.   state->x0 = s % 2048;
  193.   state->x1 = (s - state->x0) / 2048;
  194. }
  195.  
  196. static const gsl_rng_type slatec_type =
  197. {"slatec",            /* name */
  198.  4194303,            /* RAND_MAX */
  199.  0,                /* RAND_MIN */
  200.  sizeof (slatec_state_t),
  201.  &slatec_set,
  202.  &slatec_get,
  203.  &slatec_get_double};
  204.  
  205. const gsl_rng_type *gsl_rng_slatec = &slatec_type;
  206.